Stabilization of high-order solutions of the cubic Nonlinear Schrodinger Equation 
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In this paper we consider the stabilization of non- fundamental unstable stationary solutions of the 
cubic nonlinear Schrodinger equation. Specifically we study the stabilization of radially symmetric 
solutions with nodes and asymmetric complex stationary solutions. For the first ones we find partial 
stabilization similar to that recently found for vortex solutions while for the later ones stabilization 
does not seem possible. 
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I. INTRODUCTION 

Nonlinear Schrodinger Equations (NLS) are one of the 
most important models of mathematical physics arising 
in a great array of contexts [l|, [2[ as for example in semi- 
conductor electronics [3L optics in nonlinear media 
photonics [ah plasmas fj], fundamentation of quantum 
mechanics pj , dynamics of accelerators Q , mean- field 
theory of Bose- Einstein condensates [9j or in biomolecule 
dynamics [lOj , to cite a few examples. 

It is well known that in multidimensional scenarios 
there may appear concentration phenomena and collapse 
depending on the initial configuration Let us write 
the model equation in the form 



^Au + g(z)\u\ 2 u, 



(1) 



on K 2 with A = d 2 /dx 2 + d 2 /dy 2 and initial data 
uo(x,y). It is well-known that, when g{z) = go < 0, 
the solutions of Eq. (JTJ) with L 2 ( 



norm 



N(u) = \\u\\ 2 



(2) 



such that go II u II 2 is larger than a critical value N may 
undergo collapse (i.e. catastrophic formation of very 
sharp gradients and concentration phenomena for u). 
However, when <7olM| 2 < Nq there cannot be collapse 

Q. 

In the context of the analysis of the propagation of 
Kerr beams in layered optical media it was first explored 
the possibility that making the nonlinear coefficient g(z) 
to oscillate with the independent variable between values 
corresponding to the collapsing and expansion regimes 
could lead to collapse supresion [Til . 

This idea was later explored [12j | in the same context 
and exp orted to the field of Bose-Einstein condensation 
13 EH- The stabilized structure was identified as 
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a pulsating Townes soliton after some rearrangement of 
initial data [HI, [l6| that is able to propagate without es- 
sential distortions for very large distances. Some rigorous 
results concerning early time collapse (i.e. the situation 
in which collapse cannot be avoided) where found in [I?) ■ 

It has been only recently that the theoretical concept 
of stabilized solitons has been demonstrated in the labo- 
ratory in Optical experiments [HI, . 

Stabilized solitons have also been studied in three- 
dimensional scenarios [HI, H3, HH ■ Also, they have been 
considered in vector media leading to the so-called sta- 
bilized vector solitons [22l . |23| . Finally, the possibility of 
existence of stabilized solitons with different nonlineari- 
ties and dimensionalities has been studied in Ref. [24|. 

The idea of stabilized solitons and nonlincarity man- 
agement has also inspired some related mathematical 
research, with a focus on the averaged and collapse- 
preventing properties [25|, [26|, [2tJ. However, a rigorous 
theoretical mathematical description of the stabilization 
process is still missing [28j . 

The reviews [HI, H3] offer a panoramic vision of the 
field of stabilized solitons. 

In this paper we want to complement the present 
knowledge on stabilized solitons with a numerical study 
of the stabilization of more complex stationary solu- 
tions of the nonlinear Schrodinger equation. Although 
the Townes soliton can be stabilized, vortices have been 
found more difficult to stabilize [31) since the periodic 
modulation of the non-linear coefficient alone is not 
enough to achieve stabilization. 

The aim of this work is to investigate the possibility of 
stabilizing different types of stationary solutions of the 
NLS equation. Specifically, we will study the stabiliza- 
tion of higher order radially symmetric solutions in scalar 
and vector media to check if the simpler structure of these 
solutions with respect to vortices allows them to be sta- 
bilized. We will also study the stabilization of more com- 
plex asymmetric solutions of the NLSE described in Rcf. 

The paper is organized as follows: in Sec. [IT] we study 
the stabilization of higher-order radially symmetric solu- 
tions in scalar layered media and find that it is not posible 
to achieve stabilization of these structures in this simple 



2 



way. In Sec. IHII we show the stabilization of these struc- 
tures by the addition of a second stabilizing component, 
i.e. in the vector case. In Sec. IIVI we consider the pos- 
sibility of stabilizing higher order asymmetric solutions. 
Finally, in Sec. [V]we summarize our conclusions. 

II. STABILIZATION OF HIGHER-ORDER 
RADIALLY SYMMETRIC SOLUTIONS IN 
SCALAR LAYERED MEDIA 

First we start by analyzing the simple case of scalar 
fields whose propagation is governed by the NLS equation 
(fT]), with g{z) — go constant. In this case, the equation 
for stationary solutions u(z,x,y) = <f>(x, y)e lXz , is 

-A0 = -iA0 + ffo H 2 0. (3) 

It is well known that Eq. (J3]) has an infinite num- 
ber of solutions with radial symmetry each having a fi- 
nite number of nodes. We will denote these solutions by 
i?o, Ri, i?2, ■•■ labelling them by their number of nodes. 
Their norms satisfy No < Ni < A^.... Ro is the ground 
state, also called the Townes soliton, which plays an im- 
portant role in the theory since No sets a critical value for 
the norm below which collapse cannot exist. Due to the 
scaling symmetry of the nonlinear Schrodinger equation, 
these solutions exist for any value of A. 

The first two higher order stationary solutions with ra- 
dial symmetry, i.e. R\ , and R2 calculated by using a stan- 
dard shooting method, are depicted in Fig. [TJ All of our 
radially symmetric solutions to be presented in this pa- 
per have been computed by this method and then inter- 
polated onto a two-dimensional rectangular grid. After 
this injection, we have used a Newton relaxation method 
in order to increase the accuracy of the computed sta- 
tionary solutions. 

We have checked the "stationarity" of our numerically 
found solutions by propagating them through an homo- 
geneous nonlinear medium subject to the perturbation 
coming from the numerical errors (both on the initial 
data and due to the roundoff error during the evolution) . 
All our numerical simulations to be presented in this pa- 
per have been done using a second order in time and spec- 
tral in space split-step Fourier algorithm with absorbing 
boundary conditions to get rid of the outgoing radiation. 
In Fig. [2] we present the propagation of the norm and 
amplitude taking u(x,y, 0) = Ri(r). As the stationary 
solutions of Eq. ^ are all unstable in the context of Eq. 
(fT]), one expects that sooner or later the instability will 
set in. In the case of the Townes and vortex solitons the 
instability sets in by z ~ 15 and z ~ 3, respectively [3l| . 
From the propagation of the norm and amplitude shown 
in Fig. [2 we can guess that, in our case, a collapse de- 
stroying the structure of this stationary state occurs near 
z ~ 11. 

We have tried to apply the stabilization technique 
based on the modulation of the nonlinear coefficient 
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FIG. 1: Plot of Ri (solid line) and R2 (dashed line), nor- 
malized to unity, for go = —0.5, A = 0.5. For these solutions 
Ni ~ 77.17 and N 2 ~ 195.84, respectively. 
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FIG. 2: Propagation of the norm given by Eq. ((2)1 and max- 
imum amplitude A(z) = maxr X:V ) \u(x,y,z)\ for the solution 
of Eq. (TJJ with initial data u(x,y,0) = Ri(r) through an 
homogeneous nonlinear media, i.e. g(z) = —0.5. 



to the first radially excited solution by setting g(z) — 
9o + ,9i cos(r2z), with go — —0.5, <?i = —1.5 and f2 = 100 
(the parameters were set according to criteria established 
in Ref. [H [lj]) and taking u(x,y,0) = Ri(r). This 
mechanism allows a stable propagation of the Townes 
soliton over long distances of more than 400 units in adi- 
mensional units [15| but the same mechanism enhances 
only slightly the stability of singly-charged vortex soli- 
tons [3l| . In principle our stationary solution has a sim- 
pler structure than the vortex and one could expect bet- 
ter results, but we will see that this is not true. 

Our results are shown in Fig. [3J We can see how the 
inclusion of this periodic modulation leads to a shorter 
stable propagation distance of the structure and that col- 
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lapse occurs by z ~ 6.5 



described by the following set of coupled NLS equations: 




FIG. 3: Propagation of the norm given by Eq. @ and maxi- 
mum amplitude A(z) = maxjj. s ) \u(x,y, z)\ for the solution 
of Eq. (JTJ) with initial data u(x,y,0) = Ri(r) under the 
effect of a periodic modulation of the nonlinear coefficient 
g(z) = go + gi cos(f2z) where go = —0.5, g\ — —1.5 and 
Q = 100. 



We have tried to stabilize this configuration by choos- 
ing different parameters for the modulation with similar 
results. We think that the reduction of the lifetime of the 
unstable stationary structure can be understood by con- 
sidering the exchange of energy, i.e. energy flow, between 
the substructures of the first excited radially symmetric 
stationary solution: the central peak and its surround- 
ing ring. During the defocusing (g(z) > 0) stages both 
substructures spread and overlap in the region where the 
field was zero previously. Then, in the focusing (g(z) < 0) 
stage the energy from the overlapping region flows mainly 
to the central peak. Therefore, with each focusing step 
the central peak is supplied with more and more energy, 
and this process leads finally to appearance of collapse 
(see the sharp amplitude peaks for z ~ 7 in Fig. [3]). 



III. STABILIZATION OF HIGHER-ORDER 
RADIALLY SYMMETRIC SOLUTIONS IN 
VECTOR LAYERED MEDIA 



Motivation and model 



Our next idea is to try to stabilize the excited radi- 
ally symmetric stationary solutionss by using a stabilized 
Townes soliton as a guide in which the higher order so- 
lution could be stabilized as it was done in Ref. [3l| for 
vortices. Hence, we shift our attention to vector systems 



du 1 

i-g- = --Am + Si z ) (an|ui| 2 +ai 2 |u 2 | 2 ) «i, (4a) 
du 1 

i— — = --Au 2 + g(z) (a 2 i|wi| 2 + a 22 |w 2 | 2 ) u 2 . (4b) 
oz 2 

Defining a = N(u 2 )/N(ui) and Ui = Ui/Ni, i = 1,2 
Eqs. ([4]) become 

i—- = -|Aui + g{z)N 1 (an|ui| 2 + aai 2 \u 2 \ 2 ) «i(5a) 

i-r^ = -\^u 2 + g{z)Nx (a 21 \ui\ 2 + aa 22 \u 2 \ 2 ) u{£b) 

For simplicity, we will discard the tilde in what follows 
and take both components to be normalized. 

The parameter a is a measure of the strength of the 
interaction between both components, which could be 
accomplished experimentally by launching beams of dif- 
ferent energies. 



B. Case a = 

In the limit a — > 0, corresponding to the case when the 
norm of one of the components is much smaller than the 
other, Eqs. (0 become 



.dui 1 . |9 

i— — = --Ami + g(z)Niaii\ui\ m, 



.du 2 
oz 



-^Au 2 + g{z)N ia21 \ Ul \ 2 



u 2 . 



(6a) 
(6b) 



Eq. (f6"a|) is a scalar NLS equation with a modulated 
nonlinear coefficient and thus admits solutions in the 
form of stabilized Townes solitons. Eq. (|6b[) is a lin- 
ear Schrodinger equation for u 2 in a trapping potential 
generated by u\. First we look for stationary solutions 
of Eqs. ([5]) when g[z) = go of the form 



Mi(r,z) = 4> x {r) exp (-i\z) , 
u 2 {r,z) = <pp{r) cxp(-ifiz), 

i.e. solutions of the nonlinear eigenvalue problem 
1 



\<f> 



-A(p + g N\ax 



1 



-|A<£ + o A r ia 2 i|(/)| 2 



cp. 



(7) 

(8) 



(9a) 
(9b) 



Obviously, Eq. (|6ajl is equivalent to Eq. (|3|) and thus 
we get all of its stationary solutions, for instance, the 
Townes soliton. We will look for stationary solutions of 
Eq. (|9b|) with radial symmetry beyond the nodeless one 
ipo. This implies that the effective potential 



V(r)^a 21 g NM 2 , 



(10) 
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must support at least two bound states. The number 
of radially symmetric bound states supported by a two- 
dimensional potential is bounded by the inequality [34] ]: 



N2D,l=0 < 1- 



1 J R2 dr dr' rr' V(r) V(r') |ln (£) | 

2 Jdr~rV{r) 



(11) 



This inequality sets a necessary condition to have the 
required bound state (i.e. that with one node) for fixed 
go = — 27r by fixing a lower bound for the value of a,i\ 
required to obtain N2D,i=o > 2. In Table U we show the 
lower bounds and the numerical results for the number of 
bound states obtained for several values of the coefficient 
a2i. We can see that for au = 4 we can have two modes 
in the effective potential generated by Ro, the second one 
being a radially symmetric function with one node. 

a.21 N2d,i=o Bound states eigenvalues 



1 < 1.39 A ~ 0.575 

2 < 1.79 A ~ 1.925 

3 < 2.19 A ~ 3.558 

4 < 2.58 Ao ~ 5.335, Ai ~ 0.195 

TABLE I: Upper bound set by Eq. ffTTj) on the number 
of bound states in two dimensions with radial symmetry 
and without angular momentum of the potential V(r) = 
a,2igoNi \Ro | 2 . Ao and Ai are the approximate eigenvalues 
of the ground and first excited states of V(r), respectively, 
found numerically. 

We have computed the profile of the radially symmetric 
solution of Eq. (|6bj) with one node in <p\ which exists for 
tt2i = 4 by using a standard shooting method, which is 
shown in Fig. [4] 




FIG. 4: [Color online] Spatial radial profile (normalized to 
one and for A = 0.5) of the stationary solution of Eqs. © of 
the form (f> = Rq. (solid line) and the first excited state with 
radial symmetry ip = y>i (dashed line). 



cally varying nonlinear coefficient. During the propaga- 
tion the potential a 2 ig(;z)|i?o| 2 oscillates with frequency 
fl due to the periodic modulation of coefficient g(z), 
hence, the potential will change its behavior periodically 
from attractive to repulsive. In general, for the linear 
situation described by Eq. J6b|, one may achieve a non- 
dispersive propagation for small oscillations of g(z), i.e. 
without being necessary to change its sign. However, as 
we are interested to extend these results to the full nonlin- 
ear equations, we must constru ct g(z) with an alternating 
sign because analytical results [la . Il7j ] and computer sim- 
ulations [151 ] have revealed that this is a requirement to 
get stabilized solitons. 

Using results from the quantum mechanical theory of 
fast perturbations 35] we have chosen a set of param- 
eters, go = — 2-7T, gi = 8ir and fl = 100 for which the 
oscillating behavior of g(z) should maintain the profile of 
the excited solution in u%. Moreover, this set of param- 
eters leads to an stabilized Townes soliton in component 
Mi In Fig. [5] we can see that although the maximum 
value of the amplitude exhibits oscillations related to the 
periodic modulation of g(z), the norm remains constant 
which indicates that no radiation is emitted while the 
potential a,2ig{z)\u\\ 2 is switched from attractive to re- 
pulsive. 

Next we take u± — Rq, U2 = U\ as initial data for Eqs. 
This situation is described by the case of a ~ 0, 
which means that the energy injected into the medium by 
the second component U2 is much smaller in comparison 
with that of the first component u\ [3~ij . The dynamics of 
the most important parameters of the system is plotted 
in Figure [6[ 

We want to remark that the time evolution of \u\ | 2 ex- 
hibits high and low frequency oscillations corresponding 
to the periodic modulation of g(z) and to the internal 
dynamics of Eq. [6aJ respectively [la] . Therefore, the po- 
tential experienced by U2 includes this oscillation pattern. 
As the low frequency oscillations do not fulfill the require- 
ments [3l[ derived from the theory of fast perturbations, 
they lead to energy emission from the second component 
ti2 . The radiation emission is linked to the decreasing 
value of the U2 norm: as the outgoing radiation hits the 
boundary of the computational domain it is removed by 
an absorbing potential. Nevertheless, the first compo- 
nent evolves unperturbed as a stabilized Townes soliton 

0- 

The spatial intensity profile of the propagating beam 
U2, see Fig. [7J has similar shape to the initial one (see 
Fig. [1]). However, the maximum of the intensity distri- 
bution in Fig. E^b) is about three times smaller than the 
corresponding value of the initial one due to the energy 
loss during the propagation. 



C. Weak nonlinear coupling 



We have propagated numerically the initial data w 2 = 
if i according to Eq. (|6b[) alone but now with a periodi- 



Finally we turn on the nonlinear coupling in system 
(O by setting a = 0.1. The dynamics of the system pa- 
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FIG. 5: [Color online] Evolution of the (a) norm N(z) and (b) 
maximum value of the amplitude A(z) described by the Eq. 
()6b|) for parameter values go — 2ir, g\ = 87r and fl = 100. (c) 
Detailed view of the amplitude oscillations in the propagation 
range z £ [70, 72]. 
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FIG. 6: [Color online] Propagation of initial data of the form 
ui — Ro,U2 — Ui under Eq. §6§ with a = 0, with go = 
— 2n, gi = 8-7T and SI = 100. Shown are the norm N(z) and 
maximum amplitude A(z) for ui (blue) and U2 (red). 
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FIG. 7: [Color online] Pseudocolor plot of the intensity 
\u2 {x, y, z) | 2 for the same simulation as in Fig. |6]on the spatial 
region [-10, 10] x [-10, 10] for (a) z = and (b) z = 400. 



rameters is presented in Figure [5] One can see that the 
propagation of the excited state is drastically changed in 
comparison with the previous analyzed cases. Under the 
effect of the nonlinear coupling, both components mi and 
1*2 reshape their transverse spatial profile by emitting 
energy. After a short propagation distance, z — 25, the 
excited solution losses its spatial shape, acquiring a pro- 
file which resembles to the Townes soliton and which will 
be propagate along with the Townes soliton. We may say 
that both initial profiles decay or readjust (with energy 
loss) their shape, eventually leading to the formation of 
nodeless stabilized vector solitons [22]. Nevertheless, for 
smaller values of a, e.g. a = 0.05, we recover the be- 
haviour described in the previous subsection, with shape 
preservation accompanied by energy emission. 

We think that the mechanism described at the end 
of Section [TT| is responsable for the decay of component 
i*2 when the nonlinear coupling is large enough. The 
changes taking place in the component 1*2 are then cou- 
pled back to component u±, which leads to energy emis- 
sion, as seen by the decreasing norm of iti in Figure [5] 



IV. SOLUTIONS WITHOUT RADIAL 
SYMMETRY 

Equation ([3|) has more solutions beyond those hav- 
ing radial symmetry. For instance, Alfimov and cowork- 
ers [32J, using branching-off techniques from the theory 
of dynamical systems constructed solutions having non- 
trivial discrete rotational symmetries. Two examples of 
those solutions are shown in Fig. [TUJ In all of them we ob- 
serve that a large central peak is surrounded by smaller 
ones having the prescribed discrete symmetry plus an 
outer ring. 

Since these solutions are also unstable, their free prop- 
agation leads to collapse as shown in Figs. [Tl] and [12] 
Since their norms are much above the critical one both 
solutions suffer a fast instability to collapse, the constant 
behavior of the norm indicating the fact that outgoing 
radiation waves have no time to escape from this sys- 
tem before the instability develops. The sensitivity of 
those configurations to collapse manifests itself in the 



6 



s 1 

^ 0.5 




50 100 150 200 250 300 350 400 




50 100 150 200 250 300 350 400 

Z 




FIG. 10: [Color online] Pseudocolor plots of asymmetric so- 
lutions of Eq. (J3J) (a) Cs invariant solution with norm N$ ~ 
627.68 and (b) C4 invariant solution with norm N4 ~ 723.91. 
The spatial region shown is (x,y) € [—10, 10] x [—10, 10]. 



FIG. 8: [Color online] Propagation of initial data of the form 
Mi = Ro,U2 = Ui under Eq. (|6j with 01 — 0.1, with go — 
— 27r, gi = 87r and Q, — 100. Shownare the norm N(z) and 
maximum amplitude A(z) for ui (blue) and 112 (red). 
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FIG. 9: Pseudocolor plot of the intensity 112(2, y, z)\ 2 for the 
same simulation as in Fig. [5] on the spatial region for (a) 
z — on the spatial region [—4, 4] x [-4, 4] and (b) z = 400. 
on the spatial region [—2, 2] x [—2, 2]. 
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FIG. 11: [Color online] Propagation through an homogeneous 
nonlinear medium of the norm N(z) and maximum amplitude 
A(z) for the stationary solution depicted in Fig. [TOT a). 



V. CONCLUSIONS 



fact that the instability distances are two orders of mag- 
nitude smaller than those of simpler solutions such as the 
Townes, vortex or first radially excited solitons. 

We have tried to use the modulation of the nonlinear 
coefficient to stabilize the stationary solutions presented 
in Fig. |TD] However, the time evolution of the beam 
parameters does not change appreciably and in partic- 
ular, the emergence of collapse cannot be delayed. We 
have tried unsucessfully different sets of parameters to 
achieve stabilization. Intuitively, it seems difficult to be 
able to stabilize those complicated structures because of 
the coexistence of different types of substructures (peaks, 
rings, etc) and the fact that their norm ar many times 
the critical one (see the caption of Fig. [TU)) . 



In this paper we have complemented previous knowl- 
edge on stabilized solitons of the Nonlinear Schrodinger 
Equation by studying numerically the possibility of sta- 
bilizing excited stationary solutions with radial symme- 
try in both (i) scalar and (ii) vectorial layered media. 
The scalar system is unable to stabilize radially excited 
states due to the internal dynamics, i.e. energy flow, of 
this state which leads to shorter stable propagation dis- 
tances of the beam and collapse when compared with the 
propagation of the same beam through homogeneous me- 
dia. In the case of vector layered media we have shown 
how weakly coupled beams when one of the components 
is chosen to be the Townes soliton and the other an un- 
stable radially excited solution can be stabilized. In this 
situation it is found that the radially excited state ra- 
diates continuously energy while preserving an intensity 
shape similar to the initial profile. 

In both scalar and vectorial layered media, the stabi- 
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cult to achieve. This fact is due to the ireversibility of 
the internal energy flows between the wave substructures, 
e.g. rings, peaks, during propagation. Additional stabi- 
lizing mechanisms, like spatially inhomogeneous nonlin- 
earities depending on the transverse variables, might help 
in stabilizing states with nontrivial structures as it hap- 
pens in the case of non-oscillating structures [36|, [13, Ha] . 

Finally, we have tried unsuccesfully to stabilize com- 
plex asymmetric solutions of Eq. (TT|) . The complex struc- 
ture of these solutions and the fact that their power is 
many times the critical power makes stabilization prob- 
ably impossible. 
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FIG. 12: [Color online] Propagation through an homogeneous 
nonlinear medium of the norm N(z) and maximum amplitude 
A(z) for the stationary solution depicted in Fig. HOf bh 



lization of solutions bearing nontrivial structure is diffi- 
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